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We propose a numerical self-consistent method for 3D classical lattice models, which 
optimizes the variational state written as two-dimensional product of tensors. The varia- 
tional partition function is calculated by the corner transfer matrix renormalization group 
(CTMRG), which is a variant of the density matrix renormalization group (DMRG). Nu- 
merical efficiency of the method is observed via its application to the 3D Ising model. 

§1. Introduction 

Variational states represented as products of local tensols have been used widely 
for calculating the ground-state energy in one-dimensional (ID) quantum systems 
and the free energy in two-dimensional (2D) classical statistical systems. In 1941 
Kramers and Wannier calculated a lower bound of the partition function of the 2D 
Ising model using the matrix product state, ^ which is a prototype of the tensor 
product state (TPS). A remarkable point in their variational approach is that the 
calculated lower bound is a good approximation for the exact partition function 2 ) 
especially in low and high temperature regions. About 30 years later Baxter for- 
mulated a variational method, which is a natural extension of the Kramers- Wannier 
approximation, by expressing the eigenstate of the transfer matrix in the form of 
TPS. 3 )" 5 )He has also shown that the integrable model can be solved using the cor- 
ner transfer matrix generated from a local tensor with infinite degree of freedom. 
Quite recently the density matrix renormalization group (DMRG) 6 )" 8 ) has provided 
a very practical numerical procedure to perform the variational computation for the 
position-dependent TPS. 9 )" 11 ) 

A current interest about TPS is its application to higher dimensional systems. 8 )' 12 )' 13 ^ 
For example, Niggemann et al has shown that the ground state of the 2D valence- 
bond-solid (VBS) type quantum spin model is exactly expressed as the TPS, where 
the calculation of the ground-state expectation value is reduced to that of the parti- 
tion function for the corresponding 2D vertex model. 14 -*' 15 ) Hieida et al evaluated the 



*' E-mail address: nishino@phys.sci.kobe-u.ac.jp 



typeset using VTpTfX.sty <ver.l.O> 



2 T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima, Y. Akutsu, and A. Gendiar 



correlation function of the 2D deformed VBS model by combining DMRG with the 
exact variational formulation of the system. 16 ) Also for the 3D classical models, the 
authors have been generalizing the 2D TPS, in the context of the higher-dimensional 
DMRG. 13 )' 17 )' 18 ) 

In constructing the TPS, the most important step is how to represent the local 
tensors, since it often reflects on the efficiency of the variational calculation. In 
our previous works, we calculated the variational free energy of the 3D Ising model, 
employing the interaction-round-a-face (IRF) representation of the local tensor with 
16 variational parameters. 18 ) Such a variational state is often called as the IRF-type 
TPS. 12 )' 18 ) The resulting efficiency is quite good in the off-critical region, though 
the calculated transition temperature is about 1.5% higher than that obtained by the 
Monte Carlo (MC) simulations. 19 )' 20 ) However, we have only discussed the efficiency 
of the IRF-type TPS and thus it is necessary to challenge another type of TPS. 

In this paper, we construct the vertex-type TPS in 2 dimension, by introducing 
four auxiliary variables of m-states into a local tensor. The new variational state has 
2m 4 parameters, which is twice as many as that of the IRF-type TPS when m = 2. 
Numerical efficiency of this TPS is examined through its application to a 3D vertex 
model whose thermodynamic property is equivalent to the 3D Ising model. 21 ) In 
the next section, we show the construction of the 2D TPS with auxiliary variables. 
Numerical procedures to obtain the maximum of the variational partition function is 
shown in §3, and the calculated result is explained in §4. Conclusions are summarized 
in §5. 

§2. Variational Formulation 

As an example of the 3D classical model, we consider the 3D vertex model that 
has one-to-one correspondence with the simple cubic lattice Ising model. 21 ) Let us 
start from a brief description of this vertex model. Consider a simple cubic lattice of 
the size iV x x oo in the XYZ-directions, where a 2-state spin variable that takes 
either up (+) or down (— ) are sitting on each link between the nearest lattice points. 
Thus one lattice point is surrounded by 6 spins as shown in Fig. 1. This is the unit 
of the 3D vertex model, which is called a 'vertex'. We label the spins on the vertical 
links of the vertex as s (bottom) and s (top), and those on the horizontal links as a, 
a', a", and a'". 

Statistical property of the vertex model is specified by the Boltzmann weight 
w Ss ((Jo'a"(j'") assigned to the vertex. In the following, we set the weight as 

w- ss {aa'a"a"') = ]T U§ U* U* U*, t#, t#„ , (2-1) 

x=±l 

where is unity when x = y and is e f3J + \J e 2/3J — 1 otherwise; the 3D vertex model 
has the same free energy as the 3D Ising model with the Hamiltonian Ti = Jxx' , 
where x and x' denote the neighboring Ising spin variables. 21 ) 

Let us introduce two notations in order to simplify the following expressions. 
The first one is {a}, which represents the group of horizontal spins a, a', a", and a'" 
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around a vertex. For example, the vertex weight Wg S {a a 1 a" a'") is simply written as 
w Ss {a}. The second one is the matrix representation of the Boltzmann weight 



W{a] = 



w_ 



.{a} w + _{a} 
M w__{a} 



(2-2) 



where we have regarded s and s of w Ss {a}, respectively, as column and row indices 
of the 2x2 matrix. 

The system explained above can be interpreted as the infinite stack of N x N 
layers, where each layer plays the role of the layer-to-layer transfer matrix T. Using 
the matrix expression in Eq. (2.2) and writing the vertex weight at the position 
in the layer as W{a i j} 1 we can express the transfer matrix as a two-dimensionally 
connected vertices: 



T = 



(2-3) 



where 2~Z[ CT ] denotes the configuration sum for all the spins on the horizontal links, 
that are shown by black circles in Fig. 1. (We use black marks for the spin variables 
whose configuration sums are taken.) Our interest is to calculate the variational 
partition function per layer 

- (24) 

for a given TPS and further to find out the best TPS that minimizes Af^ - ]. 



T] a- 




Fig. 2. Spin variables around a variational tensor, and the variational state constructed as 2D 
product of the tensors. 
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Now we explicitly construct the trial eigenvector as the 2D tensor product, 
as is shown in Figure 2. 12 )' 1S ) Each tensor has an apex spin s = ±1 and four m-state 
auxiliary variables rj, rj', rj", and rf" shown by cubes. Let us express the elements 
of the tensor as v s {r]}, where {??} denotes the set of auxiliary variables rj, rj' , rj", and 
r]'". As we have expressed the vertex weight as a matrix in Eq. (2.2), it is useful to 
interpret v + {rj} and v ~{r]} as components of the column vector 

V{,) = ( ) . (2-5) 

Using this notation, we construct the 2D TPS as 

i*>=e n v ^ ( 2 - 6 ) 

[r?] l<ij<N 

where {f]ij} denotes the set of auxiliary variables around the lattice point and 
the configuration sum is taken over for all the auxiliary variables shown by black 
cubes in Fig. 2. Since the above construction of the TPS is similar to that of the 
transfer matrix T of the 3D vertex model, we call the TPS (2-6) as the vertex-type 
TPS. 12 )' 18 ) 




G° G 1 

Fig. 3. The vertex weight G of the 2-layer classical system (Eq. (2.9)), and G 1 of the 3-layer one 
(Eq. (2.10)). 

The main profit of expressing the variational state in the form of the TPS is that 
both and (^\T\^) can be expressed as the partition functions of 2- and 3-layer 

2D classical lattice models, respectively. 14 )' 16 )" 18 ) Let us consider (^|^) first. From 
the definition of \&) in Eq. (2.6) we obtain 

(w=e n vom ■ v fo«} = e n [e^w]. 

H [fj] ±<ij<N [rj\ [fj] l<ij<N \s=±l J 

(2-7) 

where the configuration sum is taken over for all the spins and the auxiliary variables. 
Introducing the new notation 

G'iVijhj} = V{^} • V{ % .} = J2 } (2-8) 

s=±l 
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and regarding it as a local Boltzmann weight, can be interpreted as the par- 

tition function of a 2-layer 2D vertex model 

z o^y: n G'ivijhj}- (2-9) 

lv] In] ±<ij<N 

(See Fig. 3.) In the same manner, (tp r |T|!P r ) can be expressed as the partition function 
of a 3-layer 2D vertex model: 

Z± = E II G-U^l^l^}, (2-10) 

M [<t] [*?] i<ii<^v 

which is characterized by the local Boltzmann weight: 

GH^hj} = V{^} • (w{a l3 }V{ VtJ }) = ]T } • (2-H) 

Thus once the 2m numbers of the elements (of v s {rj}) are given, the variational 
ratio A[l^] = Z 1 /^ is easily obtained via standard numerical methods for the 2D 
classical lattice models. What we have to consider next is the way to find out the 
best tensor v s {r]} which maximizes \[&]. 

§3. Maximization of the Variational Partition Function 

In order to maximize the variational partition function A[^], we consider its 
variation: 

^E^W+E^-W (3-D 

with respect to the infinitesimal change of the local tensors 

<"{'/,,} •'"{'/,,} • <»*{'/,•.,•}• (3-2) 

When the linear dimension of the system (= N) is sufficiently large, most of the 
terms in the r.h.s. of Eq. (3.1) are very close with each other. This is because the 
boundary is far away from most of the sites and hence the system can be regarded 
to be uniform in the thermodynamic limit. Thus it is sufficient to treat the variation 
of A[<^] at the center of the system: 

^ AMs ^i} foife}+ ^S} fo * fe) ' (33) 

where v s {fj c } and v s {n c } are tensors at the center c = (N/2, N/2). 

For the convenience of writing down (5 C A[^] more explicitly, we divide Z° = 
and Z 1 = {$\T\ty) into two parts as follows: 

z°= £ G^fjMx^fjM 

{Vc}{Vc} 

z x = ]T G^nMMx^ncWM. (3-4) 

{*?c } {°c } {Vc } 



6 T. Nishino, Y. Hieida, K. Okunishi, N. Maeshima, Y. Akutsu, and A. Gendiar 



The factors G and G are the "local Boltzmann weights" at the center of the 
system, and X°{fj c \ri c } and X 1 {fj c \a c \ri c } are the rest of the system. The new tensors 
X°{f] c \ri c } and X 1 {fj c \a c \rj c } play a role of "reservoirs" in the terminology of the 
DMRG, which are defined as 

x°{v c \v c } 
xl {VcK\v c } 

where the restricted product Yl(ij)^c denotes that G°{f} c \r] c } and G 1 {fj c \a c \r] c } are 
not included in the right hand sides, and the restricted sum X^]' fry an d [a]' [rj\' 
denote spin configuration sum for all the spins except for those at the center, fj c , a c , 
and r] c . Thus the division by Eq. (3.4) is equivalent to punch out the system a the 
center; this operation is similar to puncture the system in the 'puncture renormal- 
ization group' by Martm-Delgado et al. 22 ^ 

Substituting the definitions of G° and G 1 into Eq. (3.4) and introducing two 
matrices 

A Ss {f}M= ^°fekR- s 

we can express Z° and Z 1 in the binary form 

Z° = VjAV c = ]T E v*{r! c }A,s{v c \v c }v s {v c } 

^=VjBV c = E v*{Vc}B Ss {Vc\Vc}v S {Vc} (3-7) 

s {Vc } s {Vc } 

of the 2m 4 -dimensional vector V c . Substituting the above expressions into Eq. (3.3), 
we finally obtain the stationary condition (5 c A[l^] = expressed as a 2m 4 -dimensional 
generalized eigenvalue problem 

E B, s {Vc\Vc}v s {v c } = m E A *s{%\v c }v s {v c }, (3-8) 

s {Vc } s {Vc } 

which can be abbreviated as BV C = A V c . 

This is a non- linear equation for the tensors v s {fj} and v s {i]}, because the 'ma- 
trices' A and B themselves are constructed from the tensors. Therefore, equation 
(3.8) should be interpreted as a self-consistent relation for the local tensors. A way 
to find out the solution of Eq. (3.8) is to repeat numerical substitutions as follows: 

(a) Set a certain initial value to the tensor V, which has 2m 4 elements. 

(b) Numerically calculate X° and X 1 by Eqs. (3.5). This calculation can be eas- 
ily done with the help of the corner transfer matrix renormalization group 
(CTMRG). 18 )- 21 )' 23 ) 



-En g°{vm 

In}' W (ij)^c 

-Ell G^VijWijfaj}, (3-5) 

W M' M' (ij)^c 
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(c) Create A and B by Eqs. (3.6), and apply A 1 B to the tensor to obtain V' = 
A" l BV. 

(d) Create a new variational tensor V + eV where e is a small parameter. (We set 
e = 0.1.) Use V + eV as the new variational tensor and return to (b). 

(e) Stop the calculation when A \W\ does not increase any more. 

This numerical iteration works when the matrix A is regular and positive definite. 
This condition is at least satisfied in the neighborhood of the stationary point where 
X[^] takes its maximum, but is not in general for arbitrary V. Thus, the initial choice 
of the variational tensor V is relevant to the stability of the numerical calculation. 



§4. Numerical Result 



We perform a calculation for the simplest vertex-type TPS with m = 2. Note 
that each tensor has 2m 4 = 32 elements, which is twice as many as that of the 
IRF-type TPS. 18 ^ To start the numerical calculation, we set the initial tensor v s {n} 
from the vertex weight 

v'{nrfrfr{") = w +s {a a' a" a'") + aw_ s {a a' a" a'") (4-1) 

where a ~ 1 is a constant that weakly breaks the spin inversion symmetry s — > — s 
of the initial TPS; typically we set a = 1.01. The construction of the initial tensor 
explicitly uses the fact that m = 2. We then performed the numerical self-consistent 
improvement for v s {n}, and have succeeded to reach the fixed point that satisfies 
Eq. (3.8) within 1000 iterations in the whole temperature regions. 

In order to check the quality of the obtained TPS, we observe the spontaneous 
magnetization of the 3D Ising model. In the vertex representation of the 3D Ising 
model (Eq. (2.1)), the magnetization with respect to the optimized TPS is expressed 
as 

V?QV C / VjBVr 1 V?OV c 
VJAV C \VJAVJ ~ YjBY c ' 

where the new matrix O is — similar to the matrix B — defined as 

o-ssinM = J2 xl {%KK} o 8sM 

with the modified vertex weight 

o*W= E xu§uzuzuzuz,,uz„„ 

x=±l 

which represents polarization of the Ising spin at the center of the system. 

We show the calculated result in figure 4. The black triangles indicates the 
calculated spontaneous magnetization M Y (K = [3 J) by the vertex-type TPS. For 
comparison, we also show the magnetization M IRF (K) calculated by the IRF-type 
variational formulation 18 ) (white squares) and the Monte Carlo result M MC (K) 
(line) by Talapov and Blote. 19 ^ Away from the critical point, M w (K) shows good 



(4-2) 
(4-3) 

(4-4) 
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Fig. 4. Spontaneous magnetization M v (K) of the 3D Ising model calculated via TPS when m = 2. 
(See black triangles.) For comparison, magnetization obtained by the IRF-type variational 
calculation 18 ' M mF (K) (white squares) and that by a Monte Carlo simulation 19) M MC (K) 
(curve) are shown. 

agreement with M MC (ET), but as K approaches to the critical point M V (ET) deviates 
from M MC (K). The estimated critical point by the vertex-type TPS is = 0.2203, 
which is 0.6% smaller than if c MC = 0.2216544 19 )> 20 ). We can see that the result 
for the vertex-type TPS is better than that for the IRF-type variational formulation 
Kl RF = 0.2188, which is 1.5% smaller than K^ 10 . Thus we can conclude that vertex- 
type TPS is more efficient than the IRF-type TPS. 18 ) We note that the observed 
critical behavior of M v (K) in the vicinity of the critical point is the mean-field one 

M^(K) oc \J K — KY . It should be remarked that the numerical computation time 
with the vertex formulation is the same order as that of the IRF one. 

§5. Conclusions and Discussions 

We have proposed a numerical self-consistent method for 3D classical systems, 
which optimizes the 2D vertex-type TPS with 2m 4 variational parameters. The 
method is applied to the 3D vertex model, which has the same thermodynamic 
property with the 3D Ising model, and we have confirmed that the vertex-type TPS 
gives better transition temperature than the IRF-type TPS. 

If we increase the number of states m of the auxiliary variables of the vertex- 
type TPS we will be able to further tune the vertex- type TPS. It is possible to make 
a calculations of A and B in Eq. (3.6) with a short computation of a few seconds 
up to m ~ 4. We have, however, not succeeded in finding out the solution of the 
self-consistent equation (Eq. (3.8)) stably for the cases m > 2. This is because we 
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encounter a numerical problem that the matrix A becomes singular during the self- 
consistent calculation. It is our future problem to find out a systematic way to set 
up the initial tensor so that A is always regular and positive definite. 
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